/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;
using System.Diagnostics;
using dnAnalytics.LinearAlgebra.Solvers.Preconditioners;
using NUnit.Framework;
using dnAnalytics.LinearAlgebra.Solvers;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.Math;
using System.Reflection;

namespace dnAnalytics.Tests.LinearAlgebra.Solvers.Preconditioners
{
    [TestFixture]
    [Category("Managed")]
    public sealed class IncompleteLUFactorizationTest : PreconditionerTest
    {
        private static T GetMethod<T>(IncompleteLU ilu, string methodName)
        {
            var type = ilu.GetType();
            var methodInfo = type.GetMethod(methodName,
                                            BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Instance | BindingFlags.Static,
                                            null,
                                            CallingConventions.Standard,
                                            new Type[0],
                                            null);
            var obj = methodInfo.Invoke(ilu, null);
            return (T)obj;
        }

        private static Matrix GetUpperTriangle(IncompleteLU ilu)
        {
            return GetMethod<Matrix>(ilu, "UpperTriangle");
        }

        private static Matrix GetLowerTriangle(IncompleteLU ilu)
        {
            return GetMethod<Matrix>(ilu, "LowerTriangle");
        }
        
        internal override IPreConditioner CreatePreconditioner()
        {
            return new IncompleteLU();
        }

        protected override void CheckResult(IPreConditioner preconditioner, SparseMatrix matrix, Vector vector, Vector result)
        {
            Assert.AreEqual(typeof(IncompleteLU), preconditioner.GetType(), "#01");

            IncompleteLU factorization = preconditioner as IncompleteLU;
            // Compute M * result = product
            // compare vector and product. Should be equal
            Vector product = new DenseVector(result.Count);
            matrix.Multiply(result, product);

            for (int i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(vector[i],
                                                                  product[i],
                                                                  -Precision.Magnitude(s_Epsilon)),
                              "#02-" + i.ToString());
            }
        }

        [Test]
        public void CompareWithOriginalDenseMatrix()
        {
            SparseMatrix sparseMatrix = new SparseMatrix(3);
            sparseMatrix[0, 0] = -1;
            sparseMatrix[0, 1] = 5;
            sparseMatrix[0, 2] = 6;
            sparseMatrix[1, 0] = 3;
            sparseMatrix[1, 1] = -6;
            sparseMatrix[1, 2] = 1;
            sparseMatrix[2, 0] = 6;
            sparseMatrix[2, 1] = 8;
            sparseMatrix[2, 2] = 9;
            IncompleteLU ilu = new IncompleteLU();
            ilu.Initialize(sparseMatrix);
            SparseMatrix original = GetLowerTriangle(ilu).Multiply(GetUpperTriangle(ilu)) as SparseMatrix;

            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                for (int j = 0; j < sparseMatrix.Columns; j++)
                {
                    Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(sparseMatrix[i, j], 
                                                                      original[i, j],
                                                                      -Precision.Magnitude(s_Epsilon)), 
                                  "#01-" + i.ToString() + "-" + j.ToString());
                }
            }
        }
    }
}
